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A Spectral Clustering Approach to 
Lagrangian Vortex Detection 

(to appear in Physical Review E, 2016) 

Alireza Hadjighasemj^ Daniel Karraschj^ Hiroshi Teramotoj^ and George Halleij^ 

One of the ubiquitous features of real-life turbulent flows is the existence and persistence of coher¬ 
ent vortices. Here we show that such coherent vortices can be extracted as clusters of Lagrangian 
trajectories. We carry out the clustering on a weighted graph, with the weights measuring pairwise 
distances of fluid trajectories in the extended phase space of positions and time. We then extract 
coherent vortices from the graph using tools from spectral graph theory. Our method locates all 
coherent vortices in the flow simultaneously, thereby showing high potential for automated vortex 
tracking. We illustrate the performance of this technique by identifying coherent Lagrangian vortices 
in several two- and three-dimensional flows. 


I. INTRODUCTION 

It has long been recognized that even unsteady flows 
with aperiodic time dependence admit persistent pat¬ 
terns that govern the transport of passive tracers 
Generally referred to as coherent structures, these pat¬ 
terns are often vortex-type spatial features that remain 
recognizable over times exceeding typical time scales in 
the flow. Our goal here is to systematically decompose 
trajectories in such a general flow into coherent and inco¬ 
herent families, providing a conceptual simplification of 
the underlying dynamical system. 

The majority of coherent structure identification meth¬ 
ods used in fluid dynamics continues to be Eulerian (see, 
e g., IlHZl for recent examples), concerned with features 
of the instantaneous velocity field driving the flow laisi. 
The resulting Eulerian coherent structure criteria have 
been broadly used in flow structure identification, al¬ 
though none has emerged as a definitive tool of choice. 
By their focus on the velocity field, these Eulerian cri¬ 
teria inherently depend on the reference frame in which 
they are applied m- 

By contrast, Lagrangian methods identify vortical flow 
structures based on the properties of fluid particle trajec¬ 
tories [HlllIIiHIS]. Several of these methods are frame- 
invariant and hence the structures they locate (or miss) 
are the same in all frames that translate and rotate rela¬ 
tive to each other. This invariance is especially important 
for geophysical flows which are invariably defined in the 
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rotating frame of the earth. In such flows, long lived co¬ 
herent vortices may transport fluid over great distances, 
surrounded by strongly mixing background turbulence 

mini- 

Lagrangian vortex detection approaches either seek a 
coherent material boundary to the vortex, or aim to iden¬ 
tify a coherent interior of a vortex. Goherent material 
vortex boundaries are special cases of Lagrangian co¬ 
herent structures (LGSs), the most influential material 
surfaces in the flow [3]. Within this class, Lagrangian 
vortex boundaries can either be defined as outermost 
non-filamenting, closed material surfaces (elliptic LGSs 
M US). or as outermost, closed material surfaces of 
equal material rotation [HE]. Other approaches tar¬ 
get Lagrangian vortex boundaries as locations of minimal 
curvature change [18] or as curves that maximize the vol¬ 
ume to boundary size ratio throughout advection m- 

Approaches seeking the interior of Lagrangian vortices 
have mostly been probabilistic in nature. Early tech¬ 
niques relied on the diagnostic use of relative and abso¬ 
lute dispersion [2|. Later mathematical approaches offer 
a bipartition of phase space into minimally diffusive re¬ 
gions by delineating the density evolution that can be 
characterized by the Perron-Erobenius or transfer oper¬ 
ator [20H22]- Eurther diagnostic approaches have also 
been influenced by techniques for ergodic dynamical sys¬ 
tems, such as trajectory complexity and long-term aver¬ 
ages along trajectories [231125] . 

The clustering approach developed here falls in the sec¬ 
ond category, focusing on the identification of the inte¬ 
riors of coherent Lagrangian vortices. Our method is 
unconcerned with the deformation of the boundary, re¬ 
quiring only a bulk coherence for the interior of the ma¬ 
terial vortex instead. We build on techniques developed 
over the past few decades in computer science for data 
clustering [26]. While clustering methods have already 
been used in coherent structure detection in fluid flows 
[271128], here we apply spectral clustering to a graph de¬ 
scribing the spatio-temporal evolution of a fluid. This 
approach identifies coherent vortices as clusters of La¬ 
grangian trajectories remaining close over a finite-time 
interval. As we show, our proposed method detects co¬ 
herent vortices in two- and three-dimensional flows, and 
can be extended to higher dimensional problems as well. 
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Its main advantage is that it requires a relatively low 
number of Lagrangian trajectories as an input, mak¬ 
ing it suitable for the analysis of low-resolution trajec¬ 
tory data sets (see also [191 UHl (HI for methods with a 
similar capability). During the peer-review process of 
this manuscript, we were learned about the more recent 
preprint [30], which applies a similar spectral clustering 
approach to the transfer operator framework. 

Prior definitions of coherence are tied to specific ge¬ 
ometrical requirements such as convexity [m Ei], lack 
of filamentation [14], or shape coherence [TS] of the vor¬ 
tex boundary. In contrast, our approach does not pose 
any geometrical constraint on the vortex boundary, which 
helps us to identify coherent vortices that may have non- 
convex or deformable boundaries. Unlike most other La¬ 
grangian methods [HI ng Eoi ES], which rely only on 
initial and final positions of particles, our method makes 
use of intermediate particle location information (as does 
m)- Another important feature is the ability to extract 
the a priori unknown number of coherent structures from 
the trajectory data set together with their simultaneous 
detection. This is an important prerequisite for auto¬ 
matic vortex tracking in large-scale data sets (see also 
[22])- 

Our approach is based on three basic principles: 

Principle 1. [Coherence indicator] The dynamical dis¬ 
tance between two Lagrangian particles is the distance 
between their corresponding trajectories in space-time 
over a finite time interval [to,T] of interest. 

Principle 2. [Coherent structure] A coherent structure 
is a distinguished set of Lagrangian particles which have 
mutually short dynamical distances relative to the dis¬ 
tances to particles from its complement. 

This definition adopts the notion of coherence from 
spatio-temporal clustering algorithms [33] to coherence in 
fluid flows, in a fashion similar to [28] . A typical unsteady 
fluid, however, is not a union of coherent structures. 
Rather, it is composed of coherent sets and their sur¬ 
rounding incoherent background turbulence [T] [2] . Our 
third principle makes this explicit as follows. 


the transfer operator and the direct application of 
the fuzzy C-means algorithm to trajectory data sets [28] . 
We demonstrate the applicability and effectiveness of our 


method through four examples in Section IV 


II. METHOD 

The general outline of our method is as follows. To 
solve the physical we start with a discrete sam¬ 
ple of the fluid flow and generate an abstract weighted 
graph, whose nodes correspond to Lagrangian particles 
and whose edge weights are determined according to 
Next, we apply spectral clustering to this graph, which 
is particularly suited to detect clusters in the graph ac¬ 
cording to [^together with the incoherent background, 
consistently with 


A. Input: A trajectory data set 


The essential input for our algorithm is a spatio- 
temporal trajectory data set, such as particle tracks from 
a flow experiment, drifter data from the ocean, or from 
numerical integration of a differential equation. The tra¬ 
jectory data set may be sparse or spatially non-uniform 
at the initial time. Specifically, we only assume that in 
a d-dimensional configuration space, n trajectory posi¬ 
tions {x*(t)}^_^ G are available at m discrete times 
to < ti < ... < t/e < ... < tm-i = T. This informa¬ 
tion can be stored in an n x m x d-dimensional numerical 
array, with elements := x*(t/e) G 

From this trajectory data, we define the dynamical dis¬ 
tance rij between Lagrangian particles x* and x-^ as 
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Principle 3. [Coherence vs. incoherence] Coherent 
structures are surrounded by an incoherent background 
of particles. 

Our underlines the impossibility of a simple clus¬ 
tering of a general fluid flow into coherent structures. 
Instead, we formulate the following main objective. 

Problem 1. Given a fluid domain, possibly sampled dis¬ 
cretely, and a finite time interval [to,T] of interest. And 
a partition of the fluid domain into coherent structures 
surrounded by an incoherent background. 

The rest of the paper is organized as follows. Sec¬ 
tion |n| presents our method for identifying coherent vor¬ 
tices. Section [Tn] describes the relationship of our method 
with previous methods, namely the transfer operator ap¬ 
proach uniEi], its hierarchical application [34], the ap¬ 
plication of the community detection method Infomap to 


Here j-j denotes the spatial Euclidean norm, and hence 
rij approximates the L^-norm of pairwise trajectory dis¬ 
tances. Since Euclidean coordinate transformations leave 
Euclidean distances unchanged, one readily sees that the 
pairwise distances are objective^ i.e., they remain un¬ 
changed in coordinate systems rotating and translating 
relative to each other [35]. Moreover, it is noteworthy 
that the pairwise distances remain unchanged under re¬ 
finements of the spatial resolution. 


B. Similarity graph construction 

Next, we convert the spatio-temporal data set with 
the pairwise distances r^j into a similarity graph G = 
which is specified by the set of its nodes 
V = {ui,...,u^}, the set of edges E C V x V between 
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nodes, and a similarity matrix W G which asso¬ 

ciates weights Wij to the edge eij between the nodes Vi 
and Vj. 

Specifically, the nodes of G are defined as the La- 
grangian particles, i.e., Vi = xG The edges between these 
nodes have the associated weights 

~ j") (1) 

Wij = l/vij for i j, expressing pairwise similarities 
between distinct Lagrangian particles. Other definitions 
of similarity are also possible. In general, converting dis¬ 
tance to similarity can be done via any monotonically 
decreasing function, as long as the distance function Vij 
is a metric^ i.e., it satisfies for all points in the space 
the metric axioms of identity, non-negativity, symmetry, 
and triangle inequality. This is according to the intuition 
that the ordering of graph nodes from most dissimilar to 
least dissimilar should be preserved through the similar¬ 
ity conversion. 

Extending the present similarity definition 0 to the 
diagonal of W would yield infinitely large quantities. To 
regularize IT, we set the diagonal elements to a large 
constant wu = iG ^ 1, i = 1,..., n. As we shall see later, 
the actual value of K is immaterial in our algorithm. 

The entries of IT characterize the likelihood of nodes 
Vi and Vj to be in the same coherence cluster. By con¬ 
struction, IT is nonnegative {wij > 0) and symmetric 
(IT = IT^, with the superscript T referring to matrix 
transposition). 

The degree of a node Vi e V is defined as j36] 

n 

deg(Vi) :='^Wij. 
i=i 

Subsequently, the degree matrix D is defined as the di¬ 
agonal matrix with the degrees deg('r’i) on the diagonal. 
For a subset A C T of nodes, we denote its complement 
in T by A. We measure the size of A by two different 
quantities: 

1^1 := {i; Vi e A}, 

vol(^) := y]deg(vi). 

ieA 

Here, |A| measures the size of A by its number of nodes, 
while vol(A) measures the size of A by summing over the 
weights of all edges attached to nodes in A. 

C. Graph sparsification 

For large data sets, storing all entries of the similarity 
matrix IT is prohibitive. For instance, storing n = 10^ 
elements with double precision requires 8 Terabytes of 
memory, which clearly exceeds the capacity of today’s 
typical personal computers [37] . 

To address this issue, techniques have been devel¬ 
oped to sparsify IT by retaining only elements describing 
strong enough similarity. Two widely-used approaches 


are the k-nearest neighbors and the e-neighborhood ap¬ 
proaches [38]. In the former, Wij is retained if Vj (or Vi) is 
among the k nearest neighbors of Vi (or Vj)^k ^ n. In the 
latter, Wij is retained if it exceeds a specified threshold 
e. All other Wij entries are set to zero and hence require 
no storage. Other advanced spar sificat ion approaches in¬ 
clude random sampling [39|, sampling in proportion to 
edge connectivities HOI, sampling in proportion to the 
effective resistance of an edge im, and sampling using 
relative neighborhood graphs |42]-[44]. 

Here we select the e-neighborhood approach because of 
its low computational cost. For the practical determina¬ 
tion of nearest neighbors, a number of efficient packages 
are available |45l [46] . 

D. Spectral clustering 

With the notation developed so far, our original [^can 
be re-formulated as follows. 

Problem 2. [Similarity graph clustering] Given a simi¬ 
larity graph, find a partition of the set of its nodes into 
elusters such that both of the following hold: 

1. Nodes in the same cluster are similar to each other, 
which aims to maximize the within-cluster similar¬ 
ities. 

2. Nodes in a cluster are dissimilar from those located 
in other clusters or those not included in any cluster 
(incoherent background), which aims to minimize 
the between-cluster similarities. 

These two requirements for clusters implement [^and 

respectively. A particularly efficient method to identify 
clusters is spectral clustering, which we discuss below (see 
also [38] for a review). 

1. Spectral clustering and optimal graph cuts 

Given a similarity graph G = a graph eut 

is a partition of the set of nodes V into two (or possibly 
more) subsets A and B. To such a partition, we assign 
a weight eut IT (A, B) defined as the sum of the edge 
weights between two sets A and H, i.e., 

W{A,B)-.= Y. Wi^. 

ieA,jeB 

Now, consider a subset of graph nodes with very high 
within-group similarity and with weak connections to its 
complement, such as the orange set in fig.[^ A graph cut 
separating this subset from the rest of the graph (such 
as the cut indicated by the red dashed line) then yields 
a much smaller weight cut IT(A, A) than another graph 
cut through A, which would necessarily cut some of the 
strong connections within A. 

This suggests the following minimization problem, also 
known as the mineut problem^ as a solution of [^ For a 
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FIG. 1. Undirected graph partitioning. The dashed line shows 
the solution of the problem of finding a decomposition of the 
graph into two size-balanced groups with minimal number of 
edges connecting nodes from distinct groups. 


given number k of subsets, the mincut problem is to find 
a partition Ai^ ...^Ak of V which minimizes 

k 

cut{A,,...,Ak) = ^^W{Ai,Ai). ( 2 ) 

^ i=l 

For k = 2^ the mincut problem can be solved very effi¬ 
ciently (see, e.g., [U]). In practice, however, the solution 
of the mincut problem often just separates one individual 
node (the one with weakest connections) from the rest of 
the graph. One way to circumvent this problem is to pe¬ 
nalize the smallness of sets in candidate partitions. The 
most commonly applied objective functions that imple¬ 
ment this idea are the normalized cut [48], or NCut for 
short, RatioCut [49j, MinMaxCut [50] and Cheeger ratio 
cut m- Notably, not all of these graph cut objective 
functions have solutions which satisfy both conditions in 
l(cf. m for more details). 

In this paper, we use the NCut objective function, 
whose (approximate) solutions maximize the within- 
cluster similarity and minimize the between-cluster sim¬ 
ilarity: 


NCut(Ai, ...,Ak) 


1 cut(74^, Ai) 

2 vo1(74^) 


Introducing the penalizing balancing conditions, how¬ 
ever, turns the originally simple mincut problem into an 
NP hard problem [52]. Spectral clustering is a way to 
solve relaxed versions of balanced graph cut problems. 


2. Graph Laplacian 

Shi & Malik [48] showed that the solution of the Ncut 
problem can be approximated by solutions of the gener¬ 
alized eigenproblem associated with the (unnormalized) 
graph Laplacian L = D — VF, where D is the diagonal 
degree matrix of node degrees and W is the similarity 
matrix defined earlier. 


The generalized eigenvalue problem for the graph 
Laplacian is then defined as 

Lu = XDu. (3) 

We refer to its solutions as generalized eigenvectors for 
short. Generalized eigenvectors u then offer an alter¬ 
native representation of the weighted graph data. As we 
will see in the next sections, this change of representation 
enhances the cluster-properties in the data, so that clus¬ 
ters can be easily detected in the new representation. In 
particular, the simple K-means clustering algorithm has 
no difficulties to de tect the clusters in this new represen¬ 
tation (see Section [n^ regarding K-means clustering). 

It is known from Spectral Graph Theory [36] that the 
eigenvalues solving Q satisfy 0 = Ai < ... < A^. If 
the underlying graph consists of k disconnected compo¬ 
nents (clusters with zero between-cluster similarity), then 
A = 0 is a generalized eigenvalue of multiplicity k. In that 
case, the eigenspace corresponding to this eigenvalue is 
spanned by the indicator vectors of the individual con¬ 
nected components. A perturbation argument implies 
that if the between-cluster similarities remain small, then 
the eigenvectors of the first k eigenvalues remain close 
to indicator type [38]. This enables reconstructing the 
clusters from the first k eigenvectors obtained from 
The main challenge, therefore, is to extract a meaningful 
number of clusters directly from the data, as opposed to 
postulating its value beforehand. 


E. Estimating the number of clusters by 
eigenspace analysis 

For a predetermined number k, the spectral clustering 
algorithm of Shi & Malik [48] collects the k dominant 
generalized eigenvectors as cluster indicators in a matrix 
U = (i^i,... , 14 /c) G To retrieve k from the graph 

data, we adopt here the eigengap heuristic [53] by which 

k = arg min (max (gi )), (4) 

i 

where gi = A^+i — A^ for i = 1, ...,n. In other words, k 
is simply determined as the number of eigenvalues pre¬ 
ceding the largest gap in the eigenvalue sequence. The 
presence of such a gap enables us to invoke the pertur¬ 
bation argument of the previous section, and argue that 
our graph G = (V^E^W) is a perturbation of one with k 
disconnected components. 

Expression 0 determines the number of coherent clus¬ 
ters satisfying the definition given in Section |ng Ulti¬ 
mately, however, we need to partition the graph G = 
(V^E^W) into k-\-l clusters to also account for the inco¬ 
herent cluster surrounding the coherent clusters, as cod¬ 
ified in our We refer to the last, {k + l)st cluster 
arising in this process as the noise cluster or incoherent 
cluster since it includes nodes that do not belong to any 
coherent cluster. 

Spectral gap arguments were used before in the context 
of dynamical systems (see [HI Eg [Me] for examples). 
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While the number of cluster indicators (leading singular- 
and eigenvectors) in some of these works (i.e., [54l [55] ) 
coincide with the number of coherent structures, in oth¬ 
ers (i.e., [191 mi El]) the number of cluster indicators 
differs from the n umber of coherent structures (see m 
and Section III A for more details). 


Remark 1. As discussed, we identify the number of vor¬ 
tices present in a given domain by locating the largest 
gap in the eigenvalue sequence. This implies that the 
number of eigenvalues and eigenvectors to be computed 
should be greater than the maximum number of vortices 
expected to be present in a domain. In the absence of 
intuition for the maximum number of vortices, one needs 
to conduct a full matrix decomposition instead of a par¬ 
tial decomposition. The computational cost of such a 
decomposition, however, increases dramatically with re¬ 
spect to the number of eigenvalues to be computed (see 
[58] for more information). 


ALGORITHM 1. 
Input: Similarity matrix W 


^ (cf. Section IIB) 

1. Sparsify W by using the NCut algorithm (cf. Sec¬ 


tion lie ) Remove isolated nodes, i.e., nodes with de¬ 


gree zero, from G — (V, E, W). 

2. Compute the graph Laplacian L, and solve the gener¬ 
alized eigenvalue problem Lu = \Du. 


3. Identify the number k of coherent clusters as the num¬ 
ber of eigenvalues preceding the largest gap among the 
increasingly ordered eigenvalues. Select the first k gen¬ 
eralized eigenvectors ui, ...,Uk as coherent cluster indi¬ 
cators. 


4. Assemble the matrix U = Each row of 

U corresponds to a graph node (excluding the isolated 
nodes). Apply K-means to the first k eigenvectors and 
extract A: + 1 clusters. The last cluster is the incoherent 
cluster and corresponds to the mixing region filling the 
space between coherent clusters. 

Output: Clusters Ci,..., Ck+i- 


F. Retrieving clusters from matrix U by K-means 
clustering 


G. Large-scale spectral clustering 


As a last step, we employ K-means clustering to con¬ 
vert relaxed continuous spectral vectors, corresponding 
to I/’s k columns, into a discrete cluster indicator vector 
containing the cluster assignment for each node xL 
Given the spectral vectors U G and integer K, 

K-means clustering aims to determine K points in 
called centers^ so as to minimize the mean squared dis¬ 
tance from each node to its nearest center. In 1957 Stuart 
Lloyd [59] suggested a simple iterative algorithm which 
efficiently finds a local minimum for this problem. Given 
any set of K centers, the algorithm proceeds by alternat¬ 
ing between the following two steps: 


Assignment: find each node’s nearest center and as¬ 
signs it to the corresponding cluster. 

Update: recalculate cluster centers by measuring the 
mean of all nodes included in each cluster. 


These steps repeat until no node is reassigned. Readers 
not familiar with K-means can read about this algorithm 
in numerous text books, for example see [26]. Through¬ 
out the paper, we choose the number of cluster centers 
K equal to /c + I, where the last, {k + I)st cluster cor¬ 
responds to the incoherent or noise cluster discussed in 
Section HE The K-means algorithm and its probabilis¬ 
tic counterpart (fuzzy C-means) have been used before 
to extract coherent structures either directly from a tra¬ 
jectory data set [28], or indirectly from cluster indicators 
resulted from various spectral dimensionality reduction 
algorithms 


We summarize our numerical procedure in Algo¬ 
rithm [T] _ 


For large data sets, considerable time and memory is 
required to compute and store the similarity matrix W 
and the graph Laplacian L. The most commonly used ap¬ 
proach to address this issue is graph sparsification, as dis¬ 
cussed earlier in Section jll C| From the sparse similarity 
matrix W so obtained, one determines the corresponding 
Laplacian matrix L, and calls a sparse eigenvalue solver. 

Even after the sparsification of LE, however, calculat¬ 
ing the generalized eigenvectors of the graph Laplacian 
L remains challenging with 0{n^) worst-case complex¬ 
ity m- Several authors m [60] tried to alleviate the 
problem by adapting standard eigenvalue solvers to dis¬ 
tributed architecture. Other approaches are designed to 
achieve efficiency by finding numerical approximations to 
eigenfunction problems 


Here, we adopt a low-rank matrix approximation ap¬ 
proach. The main idea is to coarse-grain the similar¬ 
ity graph G = (V^E^W)^ while keeping as much in¬ 
formation as possible from the original graph and its 
weights. To this end, we construct a bipartite graph 
Gb = (yB^Ejs^Ws) from the original similarity graph 
by uniform spatial sampling of q graph nodes, called su¬ 
pernodes ^ from n graph nodes, where q n [BH [65]. A 
bipartite graph is a graph whose set of nodes Vs admits a 
partition into two disjoint sets, A and 5, such that each 
edge connects a node in A to one in 5. As a result, no 
two nodes within A and within B are connected by an 
edge. Here, we set A as the set of all n original graph 
nodes, and B as its subset of q supernodes, considered 
as independent eopies. The weights are now defined as 
before, such that the square {n q) x (n + q) similarity 
matrix Ws of the bipartite graph can be written as 


Wb 


0 

Z 0 


(5) 
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A E B 


FIG. 2. Partitioning of a bipartite graph Gb = {Vb, Eb,Wb) 
whose set of nodes Vb is divided into two disjoint sets A and B 
such that Vb = AU B. The dashed line shows the solution of 
normalized graph cut yielding a simultaneous decomposition 
of A and B. 


where Z G is a tight similarity matrix containing 

the edge weights between all nodes and supernodes, i.e., 
between A and B. Now, one can pose the Ncut problem 
to the bipartite graph whose similarity matrix enjoys a 
simple block-structure. As shown by Dhillon [ 66 ] and 
Zha et al. m. this block-structure breaks the associ¬ 
ated Ncut problem into two parts such that the domi¬ 
nant right singular vectors of the normalized g x n tight 
similarity matrix Z = D 2 ZD^ play the role of the 
generalized eigenvectors of the graph Laplacian in Sec¬ 
tion m Here, is an n x n diagonal matrix whose 
entries are column sums of Z and D 2 is a q x q diagonal 
matrix whose entries are row sums of Z (see Appendix [B] 
for more details). 

We now summarize our algorithm for large-scale tra¬ 
jectory data sets. _ 


ALGORITHM 2. 

1. Select uniformly q supernodes from n graph nodes. 

2. Gonstruct a tight similarity matrix Z G between 

all original graph nodes and the supernodes. 

3. Given Z, form Z = . Gompute the sin¬ 

gular values and vectors of Z. Select the first k right 
singular vectors ui,... ,Uk as cluster indicators for the 
original graph. 

4. Assemble the matrix U = {ui ,..., Uk). Each row of U 
corresponds to a graph node. Apply K-means to the 
first k right singular vectors and extract k 1 clus¬ 
ters. The last cluster is the incoherent cluster and cor¬ 
responds to the mixing region Filing the space between 
coherent clusters. 

Output: Glusters Ci,..., Ck+i- 


III. RELATED PREVIOUS WORK 
A. The transfer-operator approach 

In the transfer operator-based approach [2QH22] finite¬ 
time coherent sets are defined as regions in phase space 
that minimally diffuse with the surrounding phase space 
during a finite time interval. The method builds on the 
Perron-Frobenius operator or transfer operator^ which de¬ 
scribes the evolution of material densities under the flow 
map. 

In practice, the infinite-dimensional transfer operator 
needs to be approximated by a finite-dimensional ma¬ 
trix, the transition matrix P, which is most commonly 
obtained from a partition of the flow domain {Bi). and 
the flow image {Cj)- into distinct boxes, and subse¬ 
quent computation of discrete transition probabilities: 
the transition matrix entry Pij is computed as the num¬ 
ber of particles transported from Bi to Cj, normalized 
by the total number of particles released from Bi (see 
fig-|3|. This box partitioning is also referred to as Ulam’s 
method, and introduces (numerical) diffusion at the im¬ 
plementation level [ 20 ] , 

In our context, the transition matrix P can be inter¬ 
preted as the tight similarity matrix Z of a bipartite 
graph G 23 as follows: define the first set of nodes A as the 
collection of initial boxes P^, the second set of nodes B 
as the collection of final boxes Cj, and the edge weights 
as Zij = Pjj, see fig. A similar connection to spectral 
clustering and graph cuts has been worked out earlier in 
[ 68 ] . Our presentation here, however, differs from [ 68 | 
in that we interpret the graph as a bipartite graph and 
relate it to the original references [6a|67]. 

Remark 2. The size of the resulting weight matrix de¬ 
pends on the size of the Bfs and Cj’s, as well as on the 
underlying dynamics of the system. For instance, in the 
presence of chaotic dynamics, particles released at the ini¬ 
tial time can scatter in a large domain. This, in return, 
may require a large number of boxes Cj to cover the fi¬ 
nal domain, and results in a large number of columns in 
the subsequent transition matrix. In contrast, the size of 
the weight matrix of Algorithms and depends on the 
number of tracked particles. 

With this bipartite graph construction, the optimiza¬ 
tion problem which is underlying the definition of a co¬ 
herent set in the transfer-operator setting can be refor¬ 
mulated as a clustering problem. In a (bipartite) graph 
cut, such as the one shown in fig.[^ the weight of the cut 
can be interpreted as the mass leakage of one set with its 
complement. 

As discussed in Appendix minimizing the normal¬ 
ized cut for a binary cluster indicator is NP-hard. Re¬ 
laxation of the binary cluster indicator in the real value 
domain yields the eigenvector corresponding to the sec¬ 
ond smallest eigenvalue of L as an approximate cluster 
indicator [48|. However, in order to obtain a partition of 
the graph, we need to re-transform the real-valued clus¬ 
ter indicator vector of the relaxed problem into a discrete 
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FIG. 3. Interpreting transition matrix constructed from tracer advection as tight similarity matrix Z of a bipartite graph. 



FIG. 4. (a,c) Second and third largest (left) singular vectors of the normalized transition matrix for the Bickley jet flow. (b,d) 
Gorresponding coherent sets which are obtained by searching through all possible cuts ^H]- To compute the transition matrix, 
we subdivided the domain into a grid of 400 x 120 identical boxes, and released 400 particles in each box. We then advected 
particles from to = 0 to t = 40 days. 


indicator vector. The simplest way to do this is to use 
the sign of the eigenvector as a discrete cluster indicator 
function [48]. Alternatively, one can search for a split¬ 
ting point such that the resulting partition has the best 
NCut(A, A) value [48], or apply the line-search algorithm 
of [20]. Viewing the transfer operator approach [20] as 
a bipartite spectral graph partitioning EaEi] , one can 
similarly recover discrete cluster indicator vectors from 
real-valued singular vectors. 


Figure [^ shows the second largest singular vector 
of the normalized transition matrix for the Bickley jet 
model discussed in Section IIVBI We obtain the corre¬ 
sponding binary cluster indicator by searching through 
all possible Ncuts [48]. As shown in fig. 4b, the binary 
cluster indicator so obtained highlights two coherent sets 
in which coherent vortices still remain hidden. There¬ 
fore, we step in hierarchy of increasingly ordered singular 
vectors, and search for these vortices in the third singu¬ 


lar vector, shown in fig. [^ Similarly, we extract the 
corresponding discrete-valued indicator vector by exam¬ 
ining all possible Ncuts (see fig. 4d). Figure 4d reveals 
that the yellow set, as a single entity, is composed of two 
vortices. The two vortices forming the yellow set have 
overall small mass exchange with the blue set, implying 
that the objective function is minimized. The spatial 
connectedness, however, appears to be missing in this 
solution as the yellow set is composed of two distant vor¬ 
tices. In other words, a search for sets with minimal 
mass exchange, without enforcing spatial connectedness, 
can lead to a union of coherent structures as a solution 
(see also [271169] for similar observations). In this case, 
applying K-means clustering to the collection of lead¬ 
ing singular vectors may resolve the issue. However, the 
number of coherent structures that needs to be extracted, 
may not be det ectab le anymore with the eigengap heuris¬ 
tic (see Section HE), as each singular vector highlights a 













































combination of coherent structures. In fact, the number 
of coherent structures generally does not coincide with 
the number of singular/eigenvalues preceding the largest 
eigengap (see [56], pp. 1852-1853), and therefore needs 
to be guessed or to be known a priori. 

The graph cut approach, in general, does not guarantee 
that the resulting cluster will form a connected region in 
the physical space. The connectedness constraint, how¬ 
ever, can be enforced indirectly during the construction 
of the similarity graph. Our method specifically enforces 
this constraint by measuring and penalizing distances in 
the spatio-temporal domain. As a result, unlike for the 
transfer operator method, a union of coherent structures 
is not a solution for our method. In Section |IVB 


we 


will apply our Algorithrn^ to the same Bickley jet flow 
considered earlier in fig. |4] 


B. Hierarchical partitioning of the 
transfer-operator 

In the spectral clustering community, one distinguishes 
between two approaches to detect a specified number of 
clusters in a given similarity graph using the graph cut 
procedure [38l [48l |66l [67j : two-way clustering and multi¬ 
way clustering. Our methodology presented in Section [ll| 
follows (up to the introduction of the incoherent cluster) 
the multi-way clustering approach, in which k clusters 
are retrieved from the k dominant eigenvectors at once. 

In two-way clustering, the following procedure is ap¬ 
plied recursively to generate multiple clusters: (i) com¬ 
pute the top generalized eigenvector of the unnormalized 
graph Laplacian, and (ii) bisect the graph into two sub¬ 
graphs. 

In two-way clustering, the following procedure is ap¬ 
plied recursively to generate multiple clusters: (i) com¬ 
pute the top generalized eigenvector of the unnormalized 
graph Laplacian, and (ii) bisect the graph into two sub¬ 
graphs. In the transfer-operator context, this procedure 
has been put forward in m and is stopped when the ob¬ 
tained partitions no longer satisfy a pre-specified coher¬ 
ence ratio (cf. [34] for details). In the clustering analysis 
community, two-way clustering is also found to be inef¬ 
ficient due to the fact that separate eigenvalue problems 
need to be solved repeatedly [MllTOll^. 


C. Application of fuzzy clustering to a trajectory 
data set 

Recently, Froyland & Padberg-Gehle [28] proposed a 
method based on traditional fuzzy C-means clustering 
[za [73] to identify regions of phase space that remain 
compact over a finite time interval. Specifically, they 
first build a trajectory data set X G whose rows 

are vectors (X^)^=i^...^^ containing concatenated positions 
of Lagrangian particles in time. Second, they apply the 
C-means algorithm, with a prespecified number of clus¬ 
ters K and a set of K initial starting points in 


to the trajectory data set. The result is a membership 
value describing the likelihood that a trajectory belongs 
to a cluster. Thus, each trajectory carries K membership 
values, showing the degree of belonging to each of the 
K clusters. Finally, each trajectory is assigned to only 
one cluster based on the maximum membership value 
it carries. Those trajectories carrying low membership 
values for all clusters are occasionally considered to be 
non-coherent (see [28] for more details). 

Compared with the fuzzy C-mean clustering used in 
[28] . the spectral clustering technique considers the con¬ 
nectedness of the data, whereas the C-means clustering 
method considers the compactness of the data. Fuzzy 
C-means algorithm optimizes cluster compactness by as¬ 
sessing the proximity between the uncertain data points 
assigned to the cluster and the corresponding cluster cen¬ 
ter. We note that cluster centers are not true trajectories 
of a dynamical system although they are in the trajectory 
space [28]. In contrast, our spectral clustering technique 
maximizes connectedness inside clusters and disconnect¬ 
edness between clusters at the same time by measuring 
pairwise distances between trajectories. 

As opposed to centroid-based clustering algorithms 
such as K-means or C-means, where the resulting clus¬ 
ters tend to be convex sets PHI USES], spectral clustering 
can find any cluster shape, because it has no preference 
for the shape of the cluster. This is important as we 
will show in Section IIV Cl that vortices with non-convex 
shapes are the rule rather than the exception considering 
the known vortex stirring in geophysical flows m- 

Most clustering methods including centroid-based 
methods are plagued with the problem of noisy data, i.e., 
identifying good clusters amongst noise points that just 
do not belong to any cluster [77]. In some cases, even 
a few noisy points or outliers may bias the final output 
of the algorithm In our specific context, the noise 
corrsponds to the incoherent or turbulence region itself, 
where particles do not remain compact. This implies that 
the turbulence region is not residing in a hypersphere, 
and consequently cannot be captured by adding an extra 
cluster to C-means or K-means algorithms (see [77] for 
more details). 

On the other hand, the high dimensionality of the tra¬ 
jectory dataset poses a considerable challenge to K-means 
or C-means clustering approaches. First, the curse of di¬ 
mensionality can cause slow convergence for these tra¬ 
ditional algorithms, and, second, the existence of redun¬ 
dant subspaces may not allow for the identification of the 
underlying structure in the data (cf. [78] and [79], p. 10). 


Similar to many clustering methods, the K-means or 
C-means algorithms assume that the number of clusters 
K in the dataset is known beforehand which is not nec¬ 
essarily true in real-world applications. In contrast, the 
spectral clustering can detect the right number of clus¬ 
ters automatically us ing t echniques such as the eigengap 
heuristic (cf. Section HE). 

Finally, the result of K-means or C-means cluster¬ 
ing, depends on the initial guess for the cluster centers 
[Taizg, and can reach a local minimum of the objective 
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function instead of the desired global minimum [liEo]. 
Often one restarts the procedure a number of times to 
mitigate the problem. However, when the number of clus¬ 
ters K is large, the number of times to restart K-means 
or C-means to reach an optimum can be prohibitively 
high and lead to a substantial increase in runtime (cf. 

Eq]). 


IV. RESULTS 


tances over a given time interval. 

Compared with the first approach, the second ap¬ 
proach is more accurate and more memory efficient. 
However, its parallel implementation requires communi¬ 
cation between processors, which may make the com¬ 
putation prohibitively slow. For this reason, we only 
employ the second approach in our last example, the 
Arnold-Beltrami-Childress (ABC) flow, and use the first 
approach otherwise. 


We demonstrate the implementation of Algorithms 
and on four examples to detect coherent Lagrangian 
vortices. In the first example, we consider a periodi¬ 
cally forced pendulum for which we can explicitly con¬ 
firm our results using an appropriately defined Poincare 
map. Our second example is one whose temporal com¬ 
plexity is one level higher: the Bickley jet with quasi- 
periodic time dependence [HD ED- In the third exam¬ 
ple, we detect coherent Lagrangian vortices in a quasi- 
geostrophic ocean surface flow derived from satellite- 
based sea-surface height observations [83]. Our last ex¬ 
ample is a three-dimensional velocity field, the Arnold- 
Beltrami-Childress (ABC) flow, which is an exact solu¬ 
tion of Euler’s equation [84]. This is our computationally 
most demanding example, where we deploy Algorithmic 
to reduce the graph size and the associated computa¬ 
tional cost. For the rest of the examples, we use Al¬ 
gorithm IC with the e-neighborhood graph sparsification 
approach described in Section m We notice that the 
coherent structures in our first and last examples remain 
invariant in the phase space, and hence they are in prin¬ 
ciple detectable in principle by other spectral methods 
developed specifically for steady flows and maps (see e.g.. 


To implement Algorithms and [C in the forthcom¬ 
ing examples, we use a variable-order Adams-Bashforth- 
Moulton solver (ODE113 in MATLAB) to solve the dif¬ 
ferential equations. The absolute and relative tole rances 


of the ODE solver are chosen as 10 In Section IV C 


we obtain the velocity field at any given point by interpo¬ 
lating the velocity data set using bilinear interpolation. 

The dynamic distances Vij can be computed using two 
approaches that differ in terms of memory consumption, 
suitability for parallel computation and accuracy. In the 
first approach, one builds a spatio-temporal trajectory 
data set by saving trajectory positions over m interme¬ 
diate times. One then measures pairwise distances us¬ 
ing the trapezoidal rule and sparsifies them simultane¬ 
ously. This can be done effectively using the Exhaus- 
tiveSearcher model object in MATLAB or other pack¬ 
ages, such as gHUH]. This approach is memory consum¬ 
ing but highly parallelizable. 

In the second approach, one constructs the similarity 
matrix without building any spatio-temporal trajectory 
data set. To this end, one measures pairwise distances 
concurrent with the advection of particles. Specifically, 
one defines an extra output argument inside the ODE 
function which measures and cumulates the pairwise dis- 


A. The periodically forced pendulum 

Consider the periodically forced pendulum 


Xi = X2 

±2 = — sin(xi) + ecos{t). 

For 5 = 0, the system is integrable with hyperbolic 
fixed points at (0, (2m — l)7r), and elliptic fixed points 
at (0,2m7r), where m G Z. As is well known, there are 
two heteroclinic orbits connecting each successive pair of 
hyperbolic fixed points, enclosing an elliptic fixed point, 
which is in turn surrounded by periodic orbits. These 
periodic orbits appear as closed invariant curves for the 
Poincare map V := The fixed points of the flow are 
also fixed points of V. 

Kolmogorov-Arnold-Moser (KAM) theory [86j guaran¬ 
tees the survival of most closed invariant sets for V and 
0 < 5 <C 1. Increasing the perturbation strength 5 fur¬ 
ther leads to the appearance of resonance islands (HD EH] 
and to the coexistence of regular and chaotic particle tra¬ 
jectories, as one would expect in a turbulent fluid flow 
containing coherent structures. 

Figure 7b shows these surviving invariant sets (KAM 
tori and resonance islands) of the Poincare map V ob¬ 
tained for 5 = 0.4, obtained from 800 iterations of V. 
This many iterations are required to obtain continuous- 
looking boundaries of the various coherent regions. We 
would like to capture the surviving KAM regions as co¬ 
herent clusters using Algorithmic 

To construct the pairwise dynamic distances Vij and 
subsequent similarity matrix W, we advect 90,000 par¬ 
ticles, distributed initially over a uniform grid Qq of 
300 X 300 points, from to = 0 to = 800 x 27 r. The 
spatial domain ranges from —2.6 to —0.3 in xi direction 
and from —1.2 to 1.2 in ^2 direction. We output the tra¬ 
jectory data with 3600 intermediate points, evenly spaced 
in time. Moreover, we sparsify edges from the complete 
graph representing a distance greater than e = 0.45. 

Figure shows the degree of connectivity of graph 
nodes, deg(ui), as a scalar held. We refer to this scalar 
held here and in our later examples as connectivity field. 
This held looks generally smoother than other diagnostic 
fields, such as the finite-time Lyapunov exponent [891 [90| 
or finite-size Lyapunov exponent [9T1|92] fields (see flg.jC)- 
The smoothness of the connectivity held is the result of 
two averaging processes which attenuate computational 
and in-situ measurement noises. The first averaging pro¬ 
cess happens as we integrate Euclidean distances between 
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(a) (b) (c) 


FIG. 5. Comparison of three different diagnostic fields for the periodically forced pendulum. The scalar fields are constructed for 
the same integration time T = 800 x 2n. (a) Forward-time connectivity field, (b) Forward-time FTLE field, (c) Forward-time 
FSLE field. 



EIG. 6. (a) Sorted generalized eigenvalues for graph Laplacian L for the periodically forced pendulum, (b-c) The first and 
ninth generalized eigenvectors of graph Laplacian L. Isolated points resulting from the graph sparsification are shown in white. 



EIG. 7. (a) Ten clusters extracted by K-means clustering from the first nine generalized eigenvectors of graph Laplacian L for 
the periodically forced pendulum. The tenth cluster corresponds to the chaotic sea filling the space between elliptic regions, (b) 
800 iterations of the Poincare map for the periodically forced pendulum, (c) Gomputed clusters, compared with the Poincare 
map computed for the same integration time (eight hundred iterates). 
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FIG. 8. (a) The plot depicts the runtimes of Algorithm for six different resolutions for the periodically forced pendulum. 
The runtimes represent the average CPU-times for 300 processors used in parallel in these computations. The computations 
are performed on a supercomputer with 2.7GHz Intel Xeon GPUs, (b) Glustering sensitivity with respect to the sparsification 
radius. The plot shows the average within-class similarities (for nine coherent sets) relative to the e-nearest-neighbor radius 
used to sparsify the pairwise distances Vij. 




graph nodes over time. The second averaging takes place 
once we compute i.e., when summing the edge weights 
connected to a node Vi. 


Figure |6a| shows the first 20 generalized eigenvalues as 
a function of their indices. We can see that the first 
nine eigenvalues are very close to 1, while the tenth has 
an appreciable difference, creating the largest gap in the 
eigenvalue plot. This eigengap implies that the first nine 
eigenvectors are cluster indicators from which coherent 
structures should be extracted. For example, figs. 
and show the first and ninth generalized eigenvector 
of the graph Laplacian L. 

Finally, fig. [7a| shows the ten clusters extracted by the 
K-means algorithm from the first nine generalized eigen¬ 
vectors of graph Laplacian L. The tenth cluster corre¬ 
sponds to the chaotic background filling the space be¬ 
tween the coherent clusters. In fig. the extracted 
clusters are superimposed on the Poincare map, show¬ 
ing close agreement with the Lagrangian vortices of this 
example, i.e., the elliptic islands. 


Figure shows the execution times for three major 
steps of Algorithm as a function of increasing spatial 
resolution of the graph nodes. The main computational 
bottleneck, as shown in the figure, is computing the pair¬ 
wise distances and subsequently the similarity matrix W. 
For this purpose, we utilized parallel computing tech¬ 
niques with 300 CPUs, with each processor just comput¬ 
ing a few rows/columns of the sparse similarity matrix. 
Figure 8a shows the averaged CPU-times spent on each 
processor on carrying out the particle advection, sparse 
similarity matrix construction and eigen-decomposition. 


Figure shows the sensitivity of the clustering re¬ 
sults to the choice of the neighborhood radius used to 
sparsify the pairwise distances rij. In particular, the fig¬ 
ure shows how the averaged within-class similarities of 
coherent sets change with respect to the choice of neigh¬ 
borhood radius. Figure 8b suggests the existence of a 



eigenvalue index 


FIG. 9. Sorted generalized eigenvalues for the graph Lapla¬ 
cian L for the quasiperiodic Bickley jet flow. 


critical radius below which the size and shape of clus¬ 
ters can change. This critical radius simply corresponds 
to a distance where even strong edges within coherent 
sets are affected by graph sparsification. It is important 
to choose the sparsification radius such that strong edges 
will be maintained. As a rule of thumb, we set the sparsi¬ 
fication radius such that only 5%-10% of the elements in 
the similarity matrix W will be kept. To estimate such a 
radius, one can compute the pairwise distances for a sub¬ 
sample of the original graph (e.g., 40 nodes) and choose 
the sparsification radius accordingly. 


B. Quasiperiodic Bickley jet 

Next, we consider the Bickley jet, an idealized model 
of a meandering zonal jet flanked above and below by 
counter rotating vortices mi Eg. This model consists 
of a steady background flow subject to a time-dependent 
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FIG. 11. (a) Seven clusters extracted by K-means clustering from the first six generalized eigenvectors of graph Laplacian L at 
initial time, to = 0. The seventh cluster corresponds to the mixing region filling the space between the coherent clusters, (b) 
The same clusters advected passively to the final time, t = 40 day. The complete advection sequence over 40 days is illustrated 
in the online supplemental movie Ml |93j . 


perturbation. The time-dependent Hamiltonian for this 
model reads as 


= ipo{y) + 'fpi{x,y,t), 

i’oiy) = -UQLota.nh{^), 

To 

■ 3 

fn{t) exp{iknx) , 

.n=l 

where is the steady background flow and is the 
perturbation. The constants Uq and Lq are characteristic 
velocity and characteristic length scale, respectively. For 
the following analysis, we apply the set of parameters 
used in [ST] : 

Uo = 62.66 ms“^, Lq = 1770 km, kn = 
where vq = 6371 km is the mean radius of the earth. 


tpi{x,y,t) = UoLosech^{^)^ 
^0 


For fn{t) = exp(—the time-dependent part 
of the Hamiltonian consists of three Rossby waves with 
wave numbers kn traveling at speeds c^. The amplitude 
of each Rossby wave is determined by the parameters 
Sn- Specifically, the parameter values used are: ci = 
0.1446[/o, C 2 = 0.2051/0, cg = 0.461[/o, ly = 1.77 x 10^ 
ei = 0.0075, 52 = 0.15, 53 = 0.3, = 6.371 x IO^tt, 

kn = 2n7r//^. 

To construct the dynamic distances Vij and the simi¬ 
larity matrix VF, we advect 48000 particles, distributed 
initially over a uniform grid of 400 x 120 points, from 
to = 0 to t = 40 days. The spatial domain U ranges from 
0 to 20 in X direction and from —3 to 3 in ^ direction. We 
output the trajectory data with 600 intermediate points, 
evenly spaced in time. Moreover, we sparsify edges from 
the complete graph representing a distance greater than 
e = 3. 
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In fig. we show the first 20 generalized eigenvalues 
of the graph Laplacian L with respect to their indices. 
We can observe that the largest eigengap is between the 
sixth and seventh generalized eigenvalues, signaling the 
presence of six coherent clusters in the domain. Hence, 
we extract seven clusters from the first six generalized 
eigenvectors shown in figs. 10a to lOf). The last cluster, 
as described earlier in Section II E[ corresponds to the 
incoherent region filling the space between the coherent 
vortices. The observed fuzziness of the vortex boundary 
region is due to the fact that coherent and incoherent mo¬ 
tion is-on the chosen time interval-not as distinguished 
as in the forced pendulum example considered in the pre¬ 
vious section. After all, this distinction is retrieved from 
the trajectory data, as opposed to being imposed exter¬ 
nally through some threshold, for instance. Interestingly, 
this dynamic distinction is very clear in the ocean exam¬ 
ple considered in the next section, which results in very 
pronounced cluster indicators. 

Figure |lla| shows the identified clusters at the initial 
time, and fig. |llb| shows them at the final time, confirm¬ 
ing the coherence of extracted vortices over the 40-day 
period. The complete advection sequence over 40 days is 
available in the online supplemental movie MI [93] . 


C. An ocean surface data set 

Next, we apply Algorithm to a two-dimensional un¬ 
steady velocity data set obtained from AVISO satellite 
altimetry measurements [96]. The domain of the data 
set is the Agulhas leakage in the Southern Ocean, char¬ 
acterized by large coherent eddies that pinch off from the 
Agulhas current of the Indian Ocean. 

Here, we show how our coherent Lagrangian vortex 
detection principle uncovers the material eddies over in¬ 
tegration time of 168 days, ranging from to = II January 
2006 to t = 28 June 2006. The South Atlantic ocean re¬ 
gion in question is bounded by longitudes [8.5°E, I2°E] 
and latitudes [45°S, 39°S]. The region in question is cho¬ 
sen away from the coast so that particle positions will 
be available for the entire integration time. Otherwise, 
one has to discard those particles hitting obstacles or the 
coast at some intermediate times from the computation. 
We compute the pairwise accumulative distances over a 
uniform grid of 120 x 180 points using a trajectory data 
set composed of 600 evenly spaced intermediate times. 
We sparsify edges from the complete graph representing 
a distance greater than e = I. 

Figure compares the connectivity field with the 
FTLE and FSLE fields. Note that we view the connec¬ 
tivity field as a simple visualization tool from which one 
may diagnose the existence of coherent structures before 
taking the eigendecomposition step. 

In fig. |I3a[ we show the first 20 generalized eigenvalues 
of the graph Laplacian L. We can observe that the largest 
eigengap exists between the second and third general¬ 
ized eigenvalues, signaling the presence of two coherent 
clusters in the domain, which are indicated by the corre¬ 


sponding generalized eigenvectors (see figs. I3b and I3c). 


Figure [l4a| show the coherent vortices extracted from 
the first two generalized eigenvectors of graph Laplacian 
L at initial time to = II January 2006 and final time t = 
28 June 2006 respectively. In fig. |I4b[ we confirm the 
coherence of extracted vortices by advecting them to the 
final time t = 28 June 2006. 

Interestingly, the coherent cluster shown in blue con¬ 
tains isolated poin ts located far away from the cluster 
core (see fig. I4c). The presence of isolated points in 
a given cluster, however, seems to be unphysical due to 
the continuity of fluid flows. To investigate the true na¬ 
ture of these isolated points, we repeat our computation 
with a higher resolution, over a uniform grid of 300 x 300 
points, ranging from [8.5°E, I2°E ] in l ongitudes and from 
[45°S,39°S] in latitudes (see fig. I4d). The higher reso¬ 
lution computation reveals that the previously detected 
isolated points are part of a narrow Angering emanating 
from the core of the blue cluster. This is in line with the 
known vortex stirring reported by several authors (see 
eg, for example). 

Despite the strange Angering-type appearance, the 
cluster remains highly coherent over the extraction pe¬ 
riod of 168 days. The complete advection sequence over 
168 days is illustrated in the online supplemental movies 
M2 and M3 


This example underlines that a Lagrangian vortical re¬ 
gion can have an instantaneously non-convex geometry. 
It may also, over time, absorb an initial Anger-type pro¬ 
trusion and form a convex circular boundary in the end. 
This illustrates that while requiring convexity na isa, 
lack of fllamentation [T4|, or shape coherence [18] of the 
vortex boundary may yield boundaries meeting high co¬ 
herence requirements, they will not necessarily identify 
the largest set of trajectories forming a coherent cluster. 

Finally, we repeat our computation with a sparse tra¬ 
jectory data set, composed of 57 particles distributed 
non-uniformly on an unstructured grid. Here, we select 
the number of intermediate times m, and sparsifleation 
distance e similar to our earlier computation. Figure [15] 
shows the clustering result, with fig. |14a| shown in the 
background for comparison. 


D. The ABC flow 

As a last example, we consider the steady Arnold- 
Beltrami-Childress (ABC) flow [84] 

X = A sin z C cos y, 
y = B sin x -\- A cos 
z = C sin y + B cos x, 

an exact solution of Euler’s equation. We select the pa¬ 
rameter values A = B = v^, and (7 = 1. This well- 
studied set of parameter values [TH [25l EH |97] yields six 
coherent vortices. 

We construct a high resolution graph by selecting a 
uniform grid of 120 x 120 x 120 points over the spatial 
domain ranging from 0 to 27r in x, and 2 : directions. 


















14 


-39 


-40 


-41 

<D 

■o 

.|-42 

CO 

_i 

-43 


-44 


-45 



9 10 11 12 

Longitude ° 


(a) 



,4 



-39 


-40 


-41 

o 

(D 

"O 

I -42 

CO 

_i 

-43 


-44 


-45 



9 10 11 12 

Longitude ° 


(c) 


- 4.65 

- 4.7 

- 4.75 

- 4.8 

- 4.85 

- 4.9 


FIG. 12. Comparison 
integration time T = 
field. 


of three different diagnostic fields for the ocean data set. The scalar fields are constructed for the same 
168 days, (a) Forward-time connectivity field, (b) Forward-time FTLE field, (c) Forward-time FSLE 



EIG. 13. (a) Sorted generalized eigenvalues for the graph Laplacian L for the ocean data set. (b-c) The first two generalized 
eigenvectors. 



EIG. 14. (a) Coherent vortices at initial time to = H January 2006. (b) Advected image of the vortices at the final time t = 
28 June 2006. (c) Magnification of the blue cluster shown in the first panel. The figure shows some isolated points located far 
from the cluster core, (d) Corresponding cluster obtained from a higher resolution computation, revealing that the previously 
detected isolated points are part of a narrow fingering emanating from the cluster core. The complete advection sequence over 
168 days is illustrated in the online supplemental movies M2 and M3 [9^1^ . 
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FIG. 15. Coherent vor¬ 
tices are captured at 
initial time to =11 
January 2006, with a 
sparse trajectory data 
set. Figure p'4a| is shown 
in the background for 
comparison. 



Next, we subsample the phase space uniformly on a 
coarser grid by selecting q = 1000 supernodes out of 
the 120^ nodes of the original graph, and construct the 
tight similarity matrix Z G expressing similar¬ 

ity between the q supernodes and the n nodes of the 
original graph. To construct the tight similarity ma¬ 
trix Z, we measure the dynamic distances in the lifted 
system, where trajectories can flow out of the 27r cube. 
Having the similarity matrix Z in hand, we compute 
the dominant singular values and singular vectors of 
Z = ZD^^^^. The left singular eigenvectors are 

cluster indicators for the reduced graph built upon q su¬ 
pernodes, while the right singular vectors are cluster in¬ 
dicators for the original graph. 

As the last step, we retrieve seven clusters from six 
cluster indicators using the K-means algorithm. The 
last cluster, as before, shows the incoherent region fill¬ 
ing the space between the coherent clusters or vortices. 
Figure [T^ shows the six coherent clusters which are sep¬ 
arated by the incoherent cluster. The six clusters capture 
the six known coherent structures of the ABC flow iden¬ 
tified earlier in [97] . 

Due to the existence of the spatial periodic boundary 
condition, the coherent vortices are broken into pieces in 
the initial cubic domain. In figs. |16c| and |16d[ we put to¬ 
gether the pieces of six coherent vortices, and show their 
full cylindrical geometry. The colors used in figs. |16c| 
and |16d| are consistent with those in fig. |16a[ In fig. |17[ 
the clusters are superimposed on the Poincare map show¬ 
ing close agreement between the results of the two ap¬ 
proaches. 


V. CONCLUSION 

We have developed here an approach to locate coher¬ 
ent structures based on spectral graph theory. To iden¬ 
tify coherent structures, we measure the pairwise Eu¬ 
clidean distance between Lagrangian trajectories, and 
construct an undirected weighted graph describing the 
spatio-temporal evolution of fluid flows. We then iden¬ 
tify coherent vortices as clusters of Lagrangian particles 
remaining close under the flow using two different algo¬ 
rithms. In the first algorithm, we used Shi & Malik [48] 



(d) 


FIG. 16. (a) Seven clusters extracted by K-means clustering 
(/c = 7) from the first six eigenvectors of L. The first six clus¬ 
ters correspond to six coherent vortices that were identified 
earlier in El- The chaotic sea between coherent vortices is 
the seventh cluster and appears as the void between them, 
(b) The seventh cluster that appears as the chaotic sea be¬ 
tween coherent vortices, (c)-(d) 3D vortices are reconstructed 
by putting together the coherent cluster pieces. 
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due to its solid mathematical foundation and its per¬ 
formance. However, other clustering algorithms such as 
density-based clustering approaches [98] that can incor¬ 
porate the definition of noise or incoherent cluster may be 
used alternatively. Incorporating other clustering algo¬ 
rithms, and comparing their performance for the purpose 
of Lagrangian coherent vortex identification remains a vi¬ 
able future research direction. Moreover, further work is 
needed to connect graph properties with physical or me¬ 
chanical quantities characterizing the fluid motion, be¬ 
yond the heuristic and numerical arguments given in Sec¬ 
tions [l| and lEl 
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FIG. 17. Coherent vortices extracted by Algorithmic are com¬ 
pared with the Poincare map constructed for integration time 


T = 3000. 


normalized cut to identify coherent vortices whose nodes 
on graph have large internal (external) (in-)coherence. 
We demonstrate the effectiveness of the corresponding 
Algorithm to detect Lagrangian coherent vortices in 
periodic, quasiperiodic, and unsteady two-dimensional 
flows. This includes the determination of the a priori 
unknown number of present vortices in a given domain 
using the eigengap heuristic. 

In Algorithm [C we apply a recently developed graph 
sub-sampling technique [64l [65] to handle the memory 
bottleneck associated with large-scale graphs. We apply 
AlgorithmjCin our last example, the 3D steady ABC flow, 
where we succeeded to combine high sampling resolution 
with computational efficiency. 

An advantage of our approach is that it requires a 
relatively low number of Lagrangian trajectories as in¬ 
put, making it suitable for the analysis of low-resolution 
trajectory data sets (see also [T9[ |28l [29| for similar ap¬ 
proaches designed for low numbers of Lagrangian tra¬ 
jectories). Moreover, our method is taking advantage 
of trajectories’ intermediate positions, i.e., information 
that comes in most cases without additional computa¬ 
tional cost, e.g., in time resolved trajectory data sets or 
numerical integration of velocity data sets/vector fields 
(see also [28]b 

Moreover, we argue that in fluid-like flows coherence- 
related phenomena can only be conceived in the presence 
of an incoherent background, which prohibits the parti¬ 
tioning of the fluid domain into purely coherent sets or 
regions. Here, we introduced the definition of incoherent 
cluster and partitioned the fluid domain into coherent 
and incoherent clusters, an idea that appears to be miss¬ 
ing in other similar approaches UHlEilIlH]. 

Finally, we chose spectral clustering as a tool of choice 


Appendix A: Approximating Ncut 


In this section, we recall how the NCut problem can 
be solved for the case /c = 2, which partitions the graph 
into two disjoint sets. We follow closely the arguments 


Our goal is to solve the optimization problem 

min NCut(A, A). (Al) 

First, we rewrite the problem in a more convenient form. 
Given a subset A G V we define the cluster indicator 
vector / = (/i,..., Z^)"*" G with entries 



if Vi e A, 
if Vj G A. 


(A2) 


Now, Eq. (Al) can be conveniently rewritten using the 
graph Laplacian L as 


nhn/ Lf subject to / as in (A2), D/Al, / Df = vol(E). 

This is a Rayleigh quotient, and minimizing it is of com¬ 
plexity NP-hard, since we have const rain ed / to take on 
only discrete values as described in (A2). We relax the 
problem by allowing / to take arbitrary real values { 12 - 
relaxation), to obtain: 


min Lf subject to Df±l, Df = vol(E). 
feR'^ 


After substitution of g := the problem converts 

to 

min g^D~^I^LD~^I^g subject to gLD^I^X, ||^||^ = vol(E), 
geR^ 

to which the standard Rayleigh-Ritz theorem applies, 
such that its solution g is given by the second eigen¬ 
vector of ^ Re-substituting / = 
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we see that / is the second generalized eigenvector of 
Lu = \Du. 

Similarly, we can decompose the graph into k 
partitions by using cluster indicator vectors hj = 

5 • • • 5 ^n,j ) 




^/vo\(Aj) ’ 
, 0 , 


if Vi G Aj, 
otherwise, 


i = j = 


(A3) 

Then we set the matrix H € Jinxk matrix contain- 

ing those k cluster indicator vectors as columns. Observe 
that the columns in H are ortho normal to each other, 
that is H = /, and hJLhi = CMt{Ai, Ai)/\o\{Ai). So 
we can write the problem of minimizing NCut as 


Substituting (Bl) in (B2), we get 


0 

Z 0 



(l-A) 


0 




(B3) 


where is an n x n diagonal matrix whose entries are 
column sums of Z and 1^2 is an g x g diagonal matrix 
whose entries are row sums of Z. Breaking the block 
matrix form into parts, Eq. ( |B3[ ) can be rewritten as: 

q2 = (1 — A)I)igi, 

Zqi = (1 — Z)D2q2‘ 


1/2 1/2 

Let b = gi and a = D 2 g 2 , and after variable sub¬ 
stitution, we have 


min THH'^LH) 


subject to H DH = /, H as in (A3) 


Relaxing the discreteness condition and substituting T = 
we obtain the relaxed problem 

min subject to T = I. 

Again, this is the standard trace minimization prob¬ 
lem, which is solved by the matrix T composed of the 
first k eigenvectors of as columns. Re¬ 
substituting H = we see that the solution H 

consists of the first k generalized eigenvectors of Lu = 
\Du. This yields the normalized spectral clustering al¬ 
gorithm according to [48] . 


= (1 - \)h, 

= (1 - \)a. 

These equations define the SVD of the normalized matrix 
Z = . Particularly, a and b are the left 

and right singular vectors and 1 — A is the corresponding 
singular value m- 


Appendix B: Bipartite spectral graph partitioning 


In this section, we briefly recall how spectral clustering 
is applied to bipartite graphs. This specification is also 
referred to as spectral co-clustering [66j|67|, and is pre¬ 
sented here in the sub-sampling terminology introduced 
in Section |II G| It applies, however, verbatim to the bi¬ 
partite transfer-operator graph. 

Let Z G be a tight similarity matrix between 

the n graph nodes and the g supernodes. To explicitly 
capture the node-supernode relationship, we consider a 
bipartite graph Gb = {Vb, Eb^Wb) whose nodes can be 
divided into two disjoint sets A and B such that internal 
edges all have zero weights, i.e., wf^ = 0 if i;f, G A or 
vf,v^ G B. The similarity matrix of the whole bipartite 
graph Wb then reads as 


Wb = 


0 

Z 0 


(Bl) 


To partition the bipartite graph, the optimization task 
can be formalized as a generalized eigenvalue problem 
with suitable relaxation, see Appendix [A] 


LbQ = {Db - WB)q = ^DBq (B2) 


where Db is the degree matrix of Wb- 
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